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We develop path-integral Monte Carlo simulations for a parabolic two-dimensional (2D) quantum 
dot containing A'^ interacting electrons in the presence of Dresselhaus and/or Rashba spin-orbit 
couplings. Our method solves in a natural way the spin contamination problem and allows for 
numerically exact finite-temperature results at weak spin-orbit coupling. For A'' < 10 electrons, 
we present data for the addition energy, the particle density, and the total spin S in the Wigner 
molecule regime of strong Coulomb interactions. We identify magic numbers sX N = <i and N = 7 
via a peak in the addition energy. These magic numbers differ both from weak-interaction and 
classical predictions, and are stable with respect to (weak) spin-orbit couplings. 
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PACS numbers: 73.21.La, 71.70.Ej, 73.20.Qt 



I. INTRODUCTION 



Interacting few-body quantum systems have attracted 
a lot of attention over the past decades. In that respect, 
quantum dots (QDs) continue to be of fundamental in- 
terest to several fields, e.g., spintronics, nanoelectronics, 
and quantum computing. QDs are small solid-state de- 
vices, typically containing a few up to several hundred 
electrons confined in all space directions P, . They can 
be fabricated and studied using different approaches and 
materials 0, S IS S 1 ^nd control over both the charge 
and the spin degree of freedom of the confined electrons 
has been reported in experiments. Here, we address QDs 
as realized in two-dimensional (2D) semiconductor de- 
vices, e.g., in ultraclean 2D electron gases. We consider 
the case of parabolically confined individual QDs, which 
is quite appropriate in most practical cases [l|. QDs of- 
fer the possibility to tune Coulomb correlations among 
electrons via external gates. 

Spin properties of quantum dots have recently entered 
the focus of research 0, 0, 0, 0| , in particular as 
they are central both to quantum computation [l^ and 
to spintronics Here spin-orbit (SO) terms have to 

be taken into account, coupling the spin dynamics to the 
orbital motion. In general, SO coupling is a relativistic 
effect, and appears to second order in the fine structure 
constant. In most materials of interest, two main mech- 
anisms may be distinguished, namely Rashba and 
Dresselhaus SO couplings. The Rashba SO coupling 
strength an is due to the surface inversion asymmetry 
present in the confinement to a 2D electron gas, and 
therefore can be tuned by external gates The Dres- 
selhaus coupling «£) is generally not tunable but can be 
important if the host crystal has no bulk inversion sym- 
metry, e.g., in zincblende semiconductors. For simplicity, 
we only consider the linear Dresselhaus term and neglect 
various additional spin-orbit contributions, see Ref. [T^ 
for an extended discussion. These contributions can in 
principle be taken into account within our approach. Fur- 
thermore, we focus on the case of vanishing magnetic 



field. 

In this work, we investigate the behavior of a closed 
parabolic few-electron QD in the presence of spin-orbit 
couplings, containing up to iV = 9 interacting electrons. 
While the numerically exact method employed here al- 
lows to study arbitrary interactions in principle, it is 
probably most useful in the regime of intermediate-to- 
strong Coulomb interactions, where a 'Wigner molecule' 
EllialiilllElIillll is formed. In terms of the 
standard electron gas parameter Ts, the studied interac- 
tions are around Vs ~ 5 to 10. Then a clear tendency 
towards Wigner crystallization can be observed, but dif- 
ferent spatial '"shells" are not yet locked relative to each 
other as happens for very large rs [23|. In fact, while for 
r.s oo, a completely classical situation is encountered 
|25l 1231 , quantum effects still play a major role for the 
'incipient' Wigner molecule of interest here. In such a 
case, many standard calculational tools, e.g., exact di- 
agonalization p^ . the Hartree-Fock approximation, the 
fixed-node |2^ or variational Monte Carlo approach, 
or density functional theory can meet various difficul- 
ties (like artificial symmetry breakings) or require ex- 
plicit justification, see Ref. for a review. In that 
situation, finite-temperature path-integral Monte Carlo 
(PIMC) simulations represent an attractive alternative 
scheme. The case of no spin-orbit coiipling has been 
studied using PIMC in Refs. [11 IH HI ill . 

Before describing our PIMC scheme and the ensu- 
ing results, let us first discuss previous theoretical ap- 
proaches to the physics of QDs in the presence of SO 
couplings. To study their effect on single-particle ener- 
gies, one may set up perturbation theory for small SO 
couplings. Due to the linear dependence on momentum 
in Eq. (j2J) below, perturbation theory starts at second 
order and gives a quadratic decrease of single-particle 
energies with increasing Rashba coupling an (or Dressel- 
haus coupfing ) US- Single-particle energy level cross- 
ings induced by the SO coupling have been discussed in 
Refs. 10, 1^ as a function of an applied magnetic 
field. Here, we are mainly concerned with many-body ef- 
fects due to the Coulomb interaction. For two electrons. 
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exact diagonalization studies have been carried out for 
rather strong Rashba couplings and weak interactions, as 
appropriate for InSb dots 1^ ■ Energy spectra were 
examined, a jump in the magnetization of the dot as a 
function of magnetic field was found [s^ , and a favoring 
of exchange over direct interactions as a consequence of 
SO interactions was discussed [s^l- Governale [s^ has 
employed spin-density functional theory for N < 16. He 
found that a very strong SO coupling leads to new peaks 
and/or the suppression of Hund's rule peaks otherwise 
present in the addition energy spectrum. An additional 
in-plane magnetic field was argued to imply paramag- 
netic behavior. Finally, in Ref. [s^ . several approxi- 
mate schemes have been employed to study SO effects 
in weakly interacting QDs with iV « 11 to 13 electrons. 
Here we provide results for iV < 10 electrons with strong 
Coulomb interactions and SO couplings, where exact di- 
agonalization techniques may not apply anymore. One 
should also note that for « 5 to 10, spin effects are 
very important but necessitate an essentially exact treat- 
ment. For a recent comparison of QMC data to density 
functional theory in that respect, see Ref. |3^ . 

In the absence of spin-orbit couplings, PIMC simu- 
lations for QDs suffer from two well-known problems, 
namely the fermionic sign problem and the spin con- 
tamination problem p^ . (The first problem can be re- 
lieved to some extent by the multilevel blocking algo- 
rithm ^3 . Here we restrict ourselves to a simpler "brute- 
force" approach.) The spin contamination problem arises 
for an = ao = because both the total QD spin S and 
its z component Sz are good quantum numbers. Now 
Sz is in practice fixed during the simulation since there 
are no spin flip terms in the hamiltonian. At finite tem- 
perature, one then arrives at the undesirable situation 
where states with different S but the same Sz contribute 
to the simulation. This considerably complicates data 
analysis and represents a well-known problem affecting 
also other schemes, e.g., Hartree-Fock calculations. In 
the presence of spin-orbit coupling, however, neither S 
nor Sz are good quantum numbers and the full space 
of all {Sz,S} becomes accessible. With increasing SO 
couplings, we find that the sign problem worsens expo- 
nentially, restricting the applicability of our approach to 
weak SO couplings. Since in applications, SO effects are 
usually weak, however, this restriction is not too severe. 

In order to eliminate the spin contamination problem 
in the limit of zero SO coupling, we may study a few 
finite but small values for the SO couplings, and then 
extrapolate aj^^jj — > 0. This allows to reliably compute, 
for instance, the addition spectrum of the dot, where 
we find peaks (corresponding to stability islands of these 
Wigner molecules) for iV = 3 and N = 7 electrons. Fur- 
thermore, we compute the dependence of the spin state, 
(S^) = S{S + 1), as a function of particle number N, 
where S is the total spin operator. 

The structure of the paper is as follows. After intro- 
ducing the model in Section^ we derive the short-time 
propagator for interacting fcrmions in a parabolic QD 



subject to either Rashba or Dressclhaus SO coupling (or 
both of them at the same time) , and discuss the numer- 
ical scheme in some detail. Numerical results are pre- 
sented in Section IIIII and we conclude in Section IIVI 
Throughout the paper, we put h = 1. 



II. MODEL AND METHOD 

We study the A^-electron hamiltonian describing a 
closed parabolic dot in a 2D electron gas. 



H 



Pi 



^ \2m,* 

i=l ^ 



^ |r. 



N 



1=1 



SO' 



(1) 

where m* is the effective electron mass and ujq is the os- 
cillator frequency. With i = 1, . . . ,N, the vectors pi, 
denote the 2D momenta and space coordinates of all N 
electrons. The Coulomb potential among the electrons 
contains screening effects of the host material via the 
dielectric constant, k. Measuring energies (lengths) in 
units of ujo {lo = ^ / \/^^n*uJo) , a dimensionless Coulomb 
interaction parameter is given by A = / (kIqujo) ■ For 
common host materials, Iq is in the range of few up to 
hundreds of nm. The confinement energy hujo is then 
typically between 0.1 up to a few meV, which allows 
directly to determine the strength of the Coulomb in- 
teraction A. Typical A currently realized experimentally 
are between 0.5 and 5, the latter value being already 
quite close to the regime studied in our paper. For A > 1, 
Coulomb repulsion starts to dominate over the kinetic 
energy, and electrons spatially arrange on shells. This 
"Wigner molecule" regime [l^ is studied in this work. 

Let us then address the SO couplings considered here. 
Typically, two different SO couplings are of paramount 
importance in semiconductor heterostructures, namely 
the (linear) Dresselhaus and the Rashba SO coupling. We 
allow for both types and consider the spin-orbit hamilto- 
nian (of the ith particle) 

Hso = CiRiPx<y'^ - Py'^'^) + aDiPxO-"^ - PyO-^) = p- A-a, 

(2) 

where the standard Pauli matrices a ~ [a^^a"^) act in 
spin space and 



A 



-an -an 



Note that both types of SO coupling can be transformed 
into each other by a unitary transformation. Hence the 
spectra for — {a^ ^ 0) and an = (ajj 7^ 0) 
coincide (for zero magnetic field), sec also Sec. IIIII below. 

Remarkably, even on the single-particle level, there 
is no closed solution to the Schrodinger equation in a 
parabolic potential subject to SO coupling of any kind. 
We mention in passing that for a cylindrical box, the 
single-particle problem has been solved analytically in 
Ref. |30|. PIMC simulations provide a powerful tool to 
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extract numerically exact results for this few-body in- 
teracting quantum system. One starts by discretizing 
imaginary time (0 < t < f3 = l/ksT) into sufficiently 
short time intervals t = P/P, where P is the (integer) 
Trotter number. For i = r, an approximate but accu- 
rate short-time propagator can be constructed via the 
Trotter-Suzuki decomposition. To that end, we split H 
into the nonintcracting part Hq (including the SO cou- 
plings) and the remaining interaction part (~ A). As- 
suming that the spin dynamics is slow on the timescale 
r, some algebra along the lines of Ref. then gives 
the short-time single-particle propagator under Hq in the 
form 



(r2,a2|e-^^°|ri,fTi) 



27r/Q sinh(cJo''") 



(3) 



X e 



-So [r2,ri 



„ ^ — zm* (r2— ri)-A(T ^ 



where Eq = m*{aj^ + a|,). The electron is described by 
its position r and the z-component a/2 of the spin. Fur- 
thermore, 5*0 denotes the standard 2D oscillator action 



S'o[r2,ri] 



(rf -I- rf) cosh(tJo''') - 2ri • r2 
2Zq sinh(woT) 



Note that in the absence of SO couplings, the exact prop- 
agator of the harmonic oscillator is reproduced. For small 
coupling constants a_D,_R and sufficiently short r, Eq. Q 
represents a very accurate short-time approximation to 
the single-particle propagator. The spin part in the prop- 
agator can be written as 



cos(a)(5crj^ 



■ sin(a) . 



where a = {qx, ay) = 'm*A^ {r2 — ri) and a = |a|. Using 
Eq. the many-body propagator follows in the form 
o{ N X N Slater determinants. Let us then denote the 
coordinate and the spin a of the ith electron on time slice 
n (where 1 < n < P) as and tTi„, respectively. With 
respect to the time direction, we have periodic bound- 
ary conditions. Including Coulomb interactions, we then 
obtain the many-particle partition function for given dis- 
cretization T = P/P in the form 



{o-j„=±}"' n,j \n=l y 



(4) 



n,i<j 



where the N x N matrix Af^") has the matrix elements 



The last term in Eq. (0J represents the Coulomb inter- 
action between all N electrons confined to the QD. If 



there is no SO coupling, detM^"^ factorizes into a spin- 
up and a spin-down part, Sz is a constant of motion, 
and the spin contamination problem arises. We mention 
in passing that the weight in the discretized path inte- 
gral is complex- valued, and one therefore may expect 
that observables have an imaginary part. However, all 
statistical averages for physical observables must have 
zero imaginary part, and this indeed we find within the 
standard stochastic error bars. The discretized canoni- 
cal many-particle partition function then allows us to 
access equilibrium observables of interest. For concrete- 
ness, we have chosen a rather low but finite temperature, 
ksT/ujo = 0.1. Furthermore, unless stated otherwise, 
simulations were carried out for ~ and interaction 
strength A = 10, which puts us into the Wigner molecule 
regime. Note that by simply replacing an —^ao, results 
for au = follow. 

The main limitation for this type of PIMC simula- 
tion comes from the fermionic sign problem. The sign 
problem generally arises when different paths that con- 
tribute to averages carry different signs, or even complex- 
valued phases, as encountered in the case of non-zero 
SO couplings. The resulting sign cancellation when sam- 
pling fermion paths then manifests itself as a very small 
signal-to-noise ratio. For instance, as a consequence of 
exchange, these phases appear when forming Slater de- 
terminants. Here, in the presence of SO couplings, the 
sign problem occurs even for a single particle. Unfortu- 
nately, as a function of SO couplings, we find an exponen- 
tial decay of the sign, see Fig. ^ This can be rationalized 
already on the single-particle level, since the propagator 
acquires a complex phase factor in the presence of SO 
couplings. The propagator then resembles a real-time 
propagator for a single particle, with time correspond- 
ing to the SO coupling. For this problem, the exponen- 
tial severity of the sign problem is well established. In 
effect, the parameter regime where reliable simulations 
are possible is limited to small-to-intcrmediate SO cou- 
plings. We note in passing that the Rashba SO coupling 
aji — (0.4 — 1.1) X 10~"^^eUm reported for the InGaAs 
dots of Ref. are about one order of magnitude larger 
than our largest value. However, the SO coupling in In- 
GaAs is also unusally strong, and our values should ap- 
ply more directly to GaAs dots. Simultaneously, the sign 
problem becomes more severe when increasing TV or low- 
ering temperature. For the chosen parameters, the sign 
problem only allows to study QDs containing iV < 10 
electrons. The average sign is {(p) > 0.001 in all cases 
reported below. 

We then compute several observables. First, the en- 
ergy of the A^-electron dot can be obtained from 



E 



N 



dlnZjP) 



(5) 



where the derivative can be explicitly carried out using 
Eq. Knowledge of the i?jv determines the addition 
energy 



A{N) = En+i-2En + En-i. 



(6) 
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A peak in the addition energy A (A'') indicates enhanced 
stability of the A^-electron dot ("magic number") 
Note that experimentally observed addition energies are 
determined by free energy differences, while we compute 
the energy. However, for the low-temperature regime 
studied here, entropic contributions are smaller by about 
one order of magnitude, and therefore our values for 
A(A^) arc of relevance to actual experiments. Another 
quantity of interest is the total spin S, which we extract 
from the definition 

{S')=S{S+1), (7) 

where brackets denote the statistical average using 
Eq. and S is the total spin operator. Finally, spa- 
tial ordering can be monitored via the electronic particle 
density, 

n(r)=f^(<5(r-r,)), (8) 

which is normalized to / drn{r) = A^. The related charge 
densities can also be accessed experimentally, e.g., via ca- 
pacitance spectroscopy or scanning tunneling microscopy 
techniques 01. 

As we work with finite discretization r. Trotter ap- 
proximation errors have to be taken into account |39l |. 
As shown in Ref. [40| . for small r, such errors vanish 
quadratically for all observables, allowing for simple and 
efficient extrapolation schemes that completely eliminate 
this finite-T error. The scaling regime is reached for 
tloq < 0.35 when computing the energy. Furthermore, 
particle or spin densities are found to reach the scal- 
ing regime already at higher r. Results shown here have 
been carefully extrapolated down to t = with a lin- 
ear regression fit, using results from several simulations 
obtained at 1/6 < tloq < 1/3. Hence no discretization er- 
rors are present, see also Ref. and error bars denote 
just the standard stochastic Monte Carlo errors. 

Spin flip moves are an essential ingredient of the algo- 
rithm. Within the PIMC we allow for spin flips as well 
as for position moves. Single-particle moves were found 
sufficient to ensure ergodicity. The average trial step size 
for position moves was adjusted to give acceptance rates 
of the order of 30%. On the other hand, typical spin-flip 
acceptance rates were much lower (several percent) and 
strongly dependent on aji/u. The possibility to change 
Sz arises from the SO coupling and can be used to cir- 
cumvent the spin contamination problem even for the 
case of no SO coupling, namely by extrapolating finitc- 
aR results (where spins can be flipped and no spin con- 
tamination problem is present) down to — s- 0. Such a 
scheme yields the energy En as well as the expectation 
value (S^) for the total spin of the many-body system at 
finite temperature. Typically, after ^ 10^ equilibration 
passes, 1.5 X 10^ MC samples were accumulated for each 
parameter set. Our code runs at a speed of up to two 
weeks (for N — 9, P — 60) for a given parameter set per 
1.5 X 10^ samples on a standard 2 GHz Xeon processor. 



We have checked our PIMC energies for N = 2 against 
finite-temperature exact diagonalization results, includ- 
ing both interactions and spin-orbit couplings. We found 
excellent agreement, validating our approach. 



III. RESULTS 

To verify that for A = 10, wc indeed have a Wigner 
molecule let us start with the radially integrated 
particle density l|HJl shown in Fig.|21 The cylindrical sym- 
metry of the QD implies that n(r) only depends on the 
modulus r = r. The plot indicates that the sixth elec- 
tron enters the center of the dot, whereas the remaining 
five electrons arrange on an outer ring in order to min- 
imize the Coulomb repulsion. More electrons are then 
added to the outermost shell. Finally, for A^ = 9, a sec- 
ond electron enters the center. This spatial shell filling 
sequence (as opposed to orbital shell filling) is typical of 
the Wigner molecule, which forms the finite-size counter- 
part of a Wigner crystal. In fact, precisely this spatial 
filling sequence has been reported from a purely classical 
analysis (in particular, disregarding spin effects) p5ll26l|. 
At higher temperatures, the Wigner molecule melts via 
thermal fiuctuations, while for lower A, it is eventually de- 
stroyed by quantum fiuctuations. Generally, we find that 
particle densities are practically independent of the SO 
coupling strengths ur or a^i, at least for an/ulo < 0.05. 

Although we consider rather small SO couplings, the 
many-body energy En can be clearly seen to decrease as 
a function of an (here, ao = 0). This trend has been 
observed in other studies as well [s^l- In Fig.O we show 
this effect for A^ = 3 and A^ = 4 electrons. 

For urIo < 0.07, the SO couphng does not signifi- 
cantly influence addition energies for A^ < 5, see Fig. 0] 
However, there is a slight increase in the addition en- 
ergy A (6), while the peak for A^ = 7 is reduced for 
an = 0.04 (dotted curve). (In particular, A(7) = 3.00(1) 
for an 0, while A(7) = 2.95(2) at an = 0.04.) Re- 
markably, the magic numbers A^ = 3 and A = 7 en- 
countered in the Wigner molecule regime are different 
from the ones for weak interactions (where a standard 
Fermi liquid phase is present). This is indicated in the 
inset of Fig. 0] for A = 1 and = 0.04. The magic 
numbers for A = 1 can be rationalized in terms of the 
subsequent flUing of energy levels (orbitals) of a 2D har- 
monic oscillator. This predicts a peak for N = 2, where 
the lowest level is filled, and another peak at A^ = 4 
refiecting Hund's rule behavior, see also Ref. [s^. For 
strong Coulomb repulsion, the peak for A = 2 is com- 
pletely absent, while A^ = 3 now corresponds to a magic 
number. A very distinct peak in the addition energy is 
observed for N = 7. To the best of our knowledge, this 
peak is not expected for weak interactions. The picture 
of filling up spatial shells in the Wigner solid phase (dis- 
cussed above) may suggest that the filled spatial shell 
configurations N = 5 and N ~ 8 represent magic num- 
bers. However, under a classical reasoning, the addition 



5 



energy should have no pronounced peaks but exhibits a 
rather smooth and monotonic decay in A(N) [25l l26l |. 
Moreover, the classical prediction for A(A^) [2a| yields 
values one order of magnitude smaller than the A(A^) 
found here. This indicates that for A = 10, despite the 
clear onset of spatial ordering, quantum effects are still 
very important and cause the magic numbers A'^ — 3 and 
N = 7. Notably, for these N, the dot is seen to be par- 
tially spin-polarized. The magic numbers for = 3 and 
N = 7 can thus be rationalized in terms of a Hund-rule 
type behavior specific to the incipient Wigner crystallized 
regime (see the discussion below and Fig.lSJ. Our findings 
are therefore characteristic for the quantum character of 
the Wigner molecule. A purely classical "Wigner solid" 
analysis [25l l2fil | is expected to apply only for extremely 
large A. Our numerical PIMC results for En (and the 
spin) as a function of A^ and the Rashba coupling are 
summarized in Table These data were all obtained for 
A 10 and kBT/ujo = 0.1 

Next we consider the dependence of the energy on the 
two types of spin couplings. Let us take A^ = 4 electrons, 
again at A = 10, and fix [a^ + aD)^^ — 0.05. We then 
study as a function of 7 = {ur — aD)/{otR -|- au), 
where — 1 < 7 < 1 tunes the relative strength of Rashba 
versus Dresselhaus coupling. The result is depicted in 
Fig. 13 showing a symmetric curve ^£'4(7) = i5(o.o) ~ 
E[aR,an)- The Symmetric shape can be explained by the 
unitary spin rotation which transforms the Rashba- into 
the Dresselhaus term and vice versa. With increasing 
I7I, the relative weights of the Dresselhaus and Rashba 
couplings increase, leading to increasing energy gains (5i?4 
over the system without SO couplings. 

Finally, we observe that spin plays a major role for 
A = 10. Our data for S as defined in Eq. Q show a 
drastic influence of temperature. For such a strongly in- 
teracting system, energies for different S are typically 
close by, and a thermal average can then result in differ- 
ent 5* as compared to ground-state results. For instance, 
for N = 2, wc find a mixture of the singlet ground state 
5 = and the excited triplet state 5 = 2, as PIMC gives 
S{S+1) = 1.14(1) at fcsT/wo = 0.1. Due to the strong 
Coulomb interaction, there is a clear trend towards par- 
tial spin polarization as a function of A, indicated in 
Fig.lHl sec also Table|3 Moreover, the thermal average for 
(S^) obtained from PIMC simulations at A = 10 increases 
monotonically with A^, while the ground-state spin [l^ 
has rather different values and shows a nonmonotonic de- 
pendence on A^. It is worth mentioning that the peaks 
in the addition energy for A^ = 3 and N = 7 correspond 
to enhanced values of S as well. Stability of the A^ elec- 
tron dot is thus connected to a tendency towards spin 
polarization, reminiscent of Hund's rule behavior. 



For A = 4, corresponding to weaker but still sizeable in- 
teractions, the dependence of the ground-state spin is 
given in the inset. The shown nonmonotonic dependence 
on A^ reflects the standard Hund's rule physics. Com- 
paring the A = 4 and A = 10 results, we observe that 
interactions tend to further spin-polarizc the dot. The 
huge thermal effects observed in the expectation value 
(S^) also indicate that unless experiments are carried 
out at extremely low temperatures, spin blockade phe- 
nomena relying on total spin selection rules will be 
thermally washed out. Finally, we note that for the SO 
couplings studied here, spin expectation values were not 
significantly affected. 



IV. CONCLUSION 

We have investigated the behavior of up to nine elec- 
trons in a quantum dot. We took into account strong 
Coulomb correlations between the electrons, and also in- 
corporated spin-orbit couplings. Our results were ob- 
tained from path-integral quantum Monte Carlo simula- 
tions. An exponential decrease of the fermionic sign is 
found with increasing SO couplings. This sign problem 
appears even for a single electron. Nevertheless, simu- 
lations are possible for weak SO couplings, where their 
inclusion can also be used to eliminate the spin contam- 
ination problem. 

We observe peaks in the addition energy spectrum for 
A^ = 3 and N = 7, which are likely to correspond to 
the stability of partially spin-polarized configurations in- 
duced by Coulomb interactions. These peaks are neither 
expected in the weak-interaction regime nor in the classi- 
cal (deep) Wigner solid, where spin effects arc negligible. 
We hope that this prediction can soon be tested exper- 
imentally. Our results were obtained in the regime of 
weak spin-orbit couplings, since otherwise numerical in- 
stabilities associated with the sign problem occur. Given 
this restriction however, PIMC offers a powerful tool to 
analyze the effects of spin-orbit couplings in strongly in- 
teracting quantum dots. We find no dramatic effects, 
but observable downward shifts in the many-body en- 
ergy that scale quadratically in the spin-orbit couplings. 
Spin-orbit couplings also affect addition energies. 
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FIG. 2: (Color online) Charge density for = 4 to iV 
electrons at orIq = 0.02 and ao = 0, X — 10 and ksT/ 
0.1. 



TABLE I: PIMC data for En and S(S + 1) = (S^) as a func- 
tion of ctR for an = 0, A = 10, and ksT/uo — 0.1. The 
an — data are taken from extrapolations. Bracketed num- 
bers denote error estimates. 
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FIG. 1: Average sign {<j}) as a function of ur (semi- 
logarithmic scale), for A*' = 3,qd = and tljo = 0.25. The 
dotted curve is a guide to the eye only. Vertical bars denote 
standard Monte Carlo error bars (one standard deviation). 



s 

c 
LU 



17.62 
17.61 
17.6 
17.59 
31.46 

S°31.45 

1^^^31.44 



1 1 1 1 i 1 1 1 1 

: * * i , 1 

( 


' 1 ' 

(a) - 
N=3 


1 '■■■■{)-..., 
1 "■■"■■< 

1 . 1 . 1 . 1 . 


1 1 1 - 

r (b) ■ 

) 

N=4 

1 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 
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FIG. 4: (Color online) Addition energy A(A'') in units of 
ujo for A = 10. Circles denote PIMC data extrapolated to 
Or = 0. Note that A*' = 3 represents a peak, since A(3) = 
3.696(9), whereas A(2) = 3.685 and A(4) = 3.419(10). The 
dotted curve connects the corresponding PIMC results for 
Or = 0.04. Inset: Same for A = 1 and small N. The dashed 
curve is a guide to the eye only. 
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FIG. 5: Energy difference SE4 — -E{o.o) — £'(a^,Qo) for A*' = 4 
and several 7. The dotted line is a guide to the eye only. 
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FIG. 6: (Color online) S{S + 1) as defined in Eq. as 
a function of for A = 10 and ksT/uo = 0.1. Dotted or 
dashed lines are guides to the eye only. Inset; Ground-state 
(T = 0, aa/D = 0) value for 5(5'+ 1) from Ref. m for A = 10 
(squares) and for A = 4 (diamonds). 



